HerbPt1 16S Figures

require(phyloseq)
require(tidyverse)
require(phyloseq)
require(reshape2)

Loading required package: reshape2

Attaching package: ‘reshape2’

The following objects are masked from ‘package:data.table’:

    dcast, melt

The following object is masked from ‘package:tidyr’:

    smiths

require(dplyr)
require(ggplot2)

Load data

ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/DMN_ests_16S.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_rare_16S.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_hel_trans_16S.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

create alphadiversity tables

alpha_div <- estimate_richness(physeq = ps_rare, measures = c("Observed", "Shannon", "Chao1"))
#pull out metadata and concatonate with alpha diversity metrics
md<-data.frame(sample_data(ps_rare))
alpha_div_md <- rownames_to_column(alpha_div, "Barcode_ID_G") %>% full_join(md) 

Joining, by = "Barcode_ID_G"

alpha_div_md$Herbicide <- factor(alpha_div_md$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

Shannon Div plots - no significant differences among herbicide treatments at any of the three time points

ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Shannon.pdf")

Saving 7.29 x 4.51 in image

aov_t1<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)

summary(aov_t1)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.0689 0.01723   1.345  0.266
Residuals   51 0.6533 0.01281               

aov_t2<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)

summary(aov_t2)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.0449 0.01123   0.865  0.491
Residuals   49 0.6355 0.01297               

aov_t3<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)

summary(aov_t3)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.0457 0.01142   0.681  0.609
Residuals   50 0.8391 0.01678               

ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transfomred to sucessfully ordinate. At this point, the transformed data will not ordinate.

ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)

ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)
ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)

ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.05, ps_trans)
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", trymax= 300, maxit=1000)

Adonis testing of herbicide treatments by time point

ps_adonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
  print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
  vegan::adonis(physeq_dist ~ Herbicide * Time, data = md_tab)
}

ps_adonis(ps_rare_sub)

Analysis of Variance Table

Response: Distances
           Df   Sum Sq    Mean Sq F value Pr(>F)
Groups      4 0.002281 0.00057033  0.4454 0.7756
Residuals 160 0.204866 0.00128041               

Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
Herbicide        4    0.7552 0.18880  1.1017 0.02657  0.036 * 
Time             2    0.3792 0.18959  1.1063 0.01334  0.079 . 
Herbicide:Time   8    1.5797 0.19747  1.1523 0.05559  0.002 **
Residuals      150   25.7059 0.17137         0.90450          
Total          164   28.4200                 1.00000          
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ps_adonis(ps_trans_sub)

Analysis of Variance Table

Response: Distances
           Df  Sum Sq   Mean Sq F value Pr(>F)
Groups      4 0.00766 0.0019139   0.279 0.8912
Residuals 165 1.13173 0.0068590               

Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
Herbicide        4    0.3876 0.096908 0.97867 0.02300  0.558  
Time             2    0.2666 0.133309 1.34628 0.01582  0.012 *
Herbicide:Time   8    0.8525 0.106561 1.07615 0.05058  0.087 .
Residuals      155   15.3481 0.099020         0.91061         
Total          169   16.8549                  1.00000         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ps_adonis(ps_dmn)

Analysis of Variance Table

Response: Distances
           Df   Sum Sq   Mean Sq F value  Pr(>F)  
Groups      4 0.019549 0.0048873  2.8407 0.02596 *
Residuals 165 0.283879 0.0017205                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab) 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
Herbicide        4    0.2226 0.055638  2.2900 0.04823  0.001 ***
Time             2    0.1044 0.052199  2.1484 0.02263  0.005 ** 
Herbicide:Time   8    0.5214 0.065169  2.6823 0.11299  0.001 ***
Residuals      155    3.7659 0.024296         0.81615           
Total          169    4.6142                  1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ordination plots DMN

ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.0860376 
Run 1 stress 0.08603731 
... New best solution
... Procrustes: rmse 9.373782e-05  max resid 0.0005426049 
... Similar to previous best
Run 2 stress 0.08610267 
... Procrustes: rmse 0.01568934  max resid 0.1018624 
Run 3 stress 0.0867401 
Run 4 stress 0.1024122 
Run 5 stress 0.08662926 
Run 6 stress 0.1030274 
Run 7 stress 0.0861028 
... Procrustes: rmse 0.01564043  max resid 0.1017144 
Run 8 stress 0.08603769 
... Procrustes: rmse 0.0001168548  max resid 0.0006317676 
... Similar to previous best
Run 9 stress 0.09837316 
Run 10 stress 0.097449 
Run 11 stress 0.0866292 
Run 12 stress 0.08674029 
Run 13 stress 0.08673986 
Run 14 stress 0.08662905 
Run 15 stress 0.08603767 
... Procrustes: rmse 0.0008397212  max resid 0.004505676 
... Similar to previous best
Run 16 stress 0.0866296 
Run 17 stress 0.08610286 
... Procrustes: rmse 0.01562942  max resid 0.1016866 
Run 18 stress 0.08610292 
... Procrustes: rmse 0.01578954  max resid 0.1021695 
Run 19 stress 0.08674769 
Run 20 stress 0.08662955 
*** Solution reached

T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_dmn$plot + theme_classic()

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1104159 
Run 1 stress 0.1227376 
Run 2 stress 0.1095259 
... New best solution
... Procrustes: rmse 0.03288404  max resid 0.233578 
Run 3 stress 0.1122108 
Run 4 stress 0.1095259 
... New best solution
... Procrustes: rmse 0.0001772101  max resid 0.0008029145 
... Similar to previous best
Run 5 stress 0.1104151 
Run 6 stress 0.109526 
... Procrustes: rmse 0.0005058391  max resid 0.002638509 
... Similar to previous best
Run 7 stress 0.1108795 
Run 8 stress 0.1095265 
... Procrustes: rmse 0.0008754685  max resid 0.003493245 
... Similar to previous best
Run 9 stress 0.1095345 
... Procrustes: rmse 0.001373785  max resid 0.006307712 
... Similar to previous best
Run 10 stress 0.1122105 
Run 11 stress 0.1108802 
Run 12 stress 0.1095262 
... Procrustes: rmse 0.0002858933  max resid 0.001423687 
... Similar to previous best
Run 13 stress 0.1104157 
Run 14 stress 0.1104155 
Run 15 stress 0.1095258 
... New best solution
... Procrustes: rmse 0.0006211951  max resid 0.002413914 
... Similar to previous best
Run 16 stress 0.110416 
Run 17 stress 0.1095266 
... Procrustes: rmse 0.0005856117  max resid 0.002863657 
... Similar to previous best
Run 18 stress 0.1183287 
Run 19 stress 0.1095259 
... Procrustes: rmse 0.0006280806  max resid 0.00242243 
... Similar to previous best
Run 20 stress 0.1122104 
*** Solution reached

T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T2.pdf")

Saving 7.29 x 4.51 in image

ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.09696447 
Run 1 stress 0.09696557 
... Procrustes: rmse 0.0005621088  max resid 0.003280265 
... Similar to previous best
Run 2 stress 0.09692312 
... New best solution
... Procrustes: rmse 0.03971579  max resid 0.2647053 
Run 3 stress 0.09784875 
Run 4 stress 0.09784865 
Run 5 stress 0.09838606 
Run 6 stress 0.09838755 
Run 7 stress 0.09696485 
... Procrustes: rmse 0.03976122  max resid 0.2641605 
Run 8 stress 0.09701754 
... Procrustes: rmse 0.008846055  max resid 0.03955505 
Run 9 stress 0.0980118 
Run 10 stress 0.09706908 
... Procrustes: rmse 0.03980545  max resid 0.2658865 
Run 11 stress 0.1008366 
Run 12 stress 0.09706926 
... Procrustes: rmse 0.03980656  max resid 0.2661383 
Run 13 stress 0.1031707 
Run 14 stress 0.09838617 
Run 15 stress 0.09800986 
Run 16 stress 0.09777433 
Run 17 stress 0.09696528 
... Procrustes: rmse 0.03978925  max resid 0.2640928 
Run 18 stress 0.09980988 
Run 19 stress 0.09785111 
Run 20 stress 0.09696545 
... Procrustes: rmse 0.0393462  max resid 0.2628007 
Run 21 stress 0.1004913 
Run 22 stress 0.09706983 
... Procrustes: rmse 0.04026309  max resid 0.2684384 
Run 23 stress 0.09696571 
... Procrustes: rmse 0.03979867  max resid 0.2654576 
Run 24 stress 0.0969422 
... Procrustes: rmse 0.003132494  max resid 0.01735126 
Run 25 stress 0.09696556 
... Procrustes: rmse 0.03991442  max resid 0.2661971 
Run 26 stress 0.09785036 
Run 27 stress 0.09696507 
... Procrustes: rmse 0.03979892  max resid 0.2642904 
Run 28 stress 0.09850873 
Run 29 stress 0.09696463 
... Procrustes: rmse 0.00565883  max resid 0.03678046 
Run 30 stress 0.09784824 
Run 31 stress 0.09786828 
Run 32 stress 0.09696596 
... Procrustes: rmse 0.03997918  max resid 0.2659135 
Run 33 stress 0.09785936 
Run 34 stress 0.09692151 
... New best solution
... Procrustes: rmse 0.001159294  max resid 0.006295302 
... Similar to previous best
*** Solution reached

T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T3.pdf")

Saving 7.29 x 4.51 in image

Ordination plots rarefied

ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1810491 
Run 1 stress 0.1821958 
Run 2 stress 0.1815253 
... Procrustes: rmse 0.06152171  max resid 0.2535702 
Run 3 stress 0.181072 
... Procrustes: rmse 0.02192656  max resid 0.1192016 
Run 4 stress 0.180942 
... New best solution
... Procrustes: rmse 0.01250884  max resid 0.07370955 
Run 5 stress 0.182172 
Run 6 stress 0.1815091 
Run 7 stress 0.1819124 
Run 8 stress 0.1833277 
Run 9 stress 0.1823832 
Run 10 stress 0.1814063 
... Procrustes: rmse 0.02964508  max resid 0.1326549 
Run 11 stress 0.1830531 
Run 12 stress 0.1837646 
Run 13 stress 0.1814152 
... Procrustes: rmse 0.03967327  max resid 0.143494 
Run 14 stress 0.1837641 
Run 15 stress 0.1810727 
... Procrustes: rmse 0.0284734  max resid 0.1316853 
Run 16 stress 0.1814075 
... Procrustes: rmse 0.02961296  max resid 0.1335373 
Run 17 stress 0.1839169 
Run 18 stress 0.1810711 
... Procrustes: rmse 0.0284339  max resid 0.131532 
Run 19 stress 0.1812618 
... Procrustes: rmse 0.01925932  max resid 0.09503523 
Run 20 stress 0.1818106 
Run 21 stress 0.1810301 
... Procrustes: rmse 0.0257806  max resid 0.1214605 
Run 22 stress 0.1813269 
... Procrustes: rmse 0.03544089  max resid 0.1387055 
Run 23 stress 0.1815286 
Run 24 stress 0.187967 
Run 25 stress 0.1810728 
... Procrustes: rmse 0.02878186  max resid 0.1314131 
Run 26 stress 0.1835689 
Run 27 stress 0.1814025 
... Procrustes: rmse 0.02945674  max resid 0.1322751 
Run 28 stress 0.1874015 
Run 29 stress 0.1909909 
Run 30 stress 0.1814146 
... Procrustes: rmse 0.04098914  max resid 0.1452078 
Run 31 stress 0.1815991 
Run 32 stress 0.181525 
Run 33 stress 0.1814023 
... Procrustes: rmse 0.02933258  max resid 0.1318085 
Run 34 stress 0.1814154 
... Procrustes: rmse 0.04114445  max resid 0.1453592 
Run 35 stress 0.1814024 
... Procrustes: rmse 0.02945012  max resid 0.1321965 
Run 36 stress 0.1815255 
Run 37 stress 0.1821853 
Run 38 stress 0.182388 
Run 39 stress 0.1821959 
Run 40 stress 0.1828976 
Run 41 stress 0.182963 
Run 42 stress 0.1812641 
... Procrustes: rmse 0.02018832  max resid 0.09662772 
Run 43 stress 0.1808181 
... New best solution
... Procrustes: rmse 0.01560944  max resid 0.09391373 
Run 44 stress 0.1808181 
... New best solution
... Procrustes: rmse 0.0001275384  max resid 0.0007245893 
... Similar to previous best
*** Solution reached

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare$plot + theme_classic()

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1534021 
Run 1 stress 0.1501222 
... New best solution
... Procrustes: rmse 0.08230186  max resid 0.290523 
Run 2 stress 0.1495892 
... New best solution
... Procrustes: rmse 0.05350526  max resid 0.2392572 
Run 3 stress 0.1502109 
Run 4 stress 0.149862 
... Procrustes: rmse 0.02211  max resid 0.07569195 
Run 5 stress 0.1531695 
Run 6 stress 0.1492979 
... New best solution
... Procrustes: rmse 0.02816165  max resid 0.1578242 
Run 7 stress 0.1547779 
Run 8 stress 0.1509426 
Run 9 stress 0.1493867 
... Procrustes: rmse 0.0124305  max resid 0.06086262 
Run 10 stress 0.1493291 
... Procrustes: rmse 0.02028147  max resid 0.1029242 
Run 11 stress 0.1493185 
... Procrustes: rmse 0.01399418  max resid 0.08297848 
Run 12 stress 0.1509814 
Run 13 stress 0.1502651 
Run 14 stress 0.1494113 
... Procrustes: rmse 0.01538676  max resid 0.07614067 
Run 15 stress 0.1508345 
Run 16 stress 0.1532662 
Run 17 stress 0.1536907 
Run 18 stress 0.1506317 
Run 19 stress 0.1533629 
Run 20 stress 0.1493236 
... Procrustes: rmse 0.01461814  max resid 0.08672645 
Run 21 stress 0.1499422 
Run 22 stress 0.1558607 
Run 23 stress 0.1508846 
Run 24 stress 0.151118 
Run 25 stress 0.152464 
Run 26 stress 0.154282 
Run 27 stress 0.1571266 
Run 28 stress 0.1500629 
Run 29 stress 0.1496188 
... Procrustes: rmse 0.01580802  max resid 0.06869795 
Run 30 stress 0.1499622 
Run 31 stress 0.1529377 
Run 32 stress 0.1543548 
Run 33 stress 0.1497197 
... Procrustes: rmse 0.02742456  max resid 0.1268267 
Run 34 stress 0.1494947 
... Procrustes: rmse 0.01488277  max resid 0.05659169 
Run 35 stress 0.1497881 
... Procrustes: rmse 0.01953659  max resid 0.07442972 
Run 36 stress 0.1496531 
... Procrustes: rmse 0.04709876  max resid 0.1934666 
Run 37 stress 0.149388 
... Procrustes: rmse 0.01158883  max resid 0.06096413 
Run 38 stress 0.1493971 
... Procrustes: rmse 0.01766645  max resid 0.1038983 
Run 39 stress 0.1502748 
Run 40 stress 0.1496127 
... Procrustes: rmse 0.02905464  max resid 0.1614178 
Run 41 stress 0.1496629 
... Procrustes: rmse 0.01445383  max resid 0.05956121 
Run 42 stress 0.1493002 
... Procrustes: rmse 0.002545233  max resid 0.01386556 
Run 43 stress 0.1493066 
... Procrustes: rmse 0.0110654  max resid 0.06490723 
Run 44 stress 0.1504065 
Run 45 stress 0.1502779 
Run 46 stress 0.1492988 
... Procrustes: rmse 0.008352935  max resid 0.04909222 
Run 47 stress 0.151426 
Run 48 stress 0.1493351 
... Procrustes: rmse 0.02178304  max resid 0.1112067 
Run 49 stress 0.1497579 
... Procrustes: rmse 0.02964951  max resid 0.1628834 
Run 50 stress 0.1493165 
... Procrustes: rmse 0.01361412  max resid 0.08043151 
Run 51 stress 0.1501407 
Run 52 stress 0.1519415 
Run 53 stress 0.1511297 
Run 54 stress 0.1497795 
... Procrustes: rmse 0.02848309  max resid 0.1510366 
Run 55 stress 0.1493392 
... Procrustes: rmse 0.02236268  max resid 0.1148816 
Run 56 stress 0.1502676 
Run 57 stress 0.1514846 
Run 58 stress 0.1497267 
... Procrustes: rmse 0.01883783  max resid 0.08695915 
Run 59 stress 0.1492971 
... New best solution
... Procrustes: rmse 0.002159199  max resid 0.009525525 
... Similar to previous best
*** Solution reached

T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2.pdf")

Saving 7.29 x 4.51 in image

ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1249587 
Run 1 stress 0.1259781 
Run 2 stress 0.1250301 
... Procrustes: rmse 0.01325505  max resid 0.06684092 
Run 3 stress 0.1251203 
... Procrustes: rmse 0.02363859  max resid 0.07522807 
Run 4 stress 0.1256432 
Run 5 stress 0.1251302 
... Procrustes: rmse 0.01336682  max resid 0.06533419 
Run 6 stress 0.1259261 
Run 7 stress 0.1259326 
Run 8 stress 0.125684 
Run 9 stress 0.1266757 
Run 10 stress 0.1257443 
Run 11 stress 0.1275958 
Run 12 stress 0.1266432 
Run 13 stress 0.1249629 
... Procrustes: rmse 0.003407392  max resid 0.01309949 
Run 14 stress 0.1250876 
... Procrustes: rmse 0.02158747  max resid 0.07005849 
Run 15 stress 0.1256465 
Run 16 stress 0.1257186 
Run 17 stress 0.1256958 
Run 18 stress 0.1250745 
... Procrustes: rmse 0.008385728  max resid 0.0374647 
Run 19 stress 0.125072 
... Procrustes: rmse 0.01857748  max resid 0.07709138 
Run 20 stress 0.1250324 
... Procrustes: rmse 0.006415311  max resid 0.03741141 
Run 21 stress 0.1257058 
Run 22 stress 0.1259713 
Run 23 stress 0.1266542 
Run 24 stress 0.1250275 
... Procrustes: rmse 0.006788055  max resid 0.03673037 
Run 25 stress 0.1250087 
... Procrustes: rmse 0.005547921  max resid 0.01757564 
Run 26 stress 0.1250586 
... Procrustes: rmse 0.01199068  max resid 0.03664954 
Run 27 stress 0.1254315 
... Procrustes: rmse 0.01544117  max resid 0.04724371 
Run 28 stress 0.1249508 
... New best solution
... Procrustes: rmse 0.001496567  max resid 0.004758168 
... Similar to previous best
*** Solution reached

T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3.pdf")

Saving 7.29 x 4.51 in image

box and whisker plots of distance within group distances

library(micrUBIfuns)
beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")

$data

$plot

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")

$data

$plot

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide")

$data

$plot

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image

treatment to control

plotDistances = function(p, m, s, d) {

  # calc distances
  wu = phyloseq::distance(p, m)
  wu.m = melt(as.matrix(wu))
  
  # remove self-comparisons
  wu.m = wu.m %>%
    filter(as.character(Var1) != as.character(Var2)) %>%
    mutate_if(is.factor,as.character)
  
  # get sample data (S4 error OK and expected)
  sd = data.frame(sample_data(p)) %>%
    select(s, d) %>%
    mutate_if(is.factor,as.character)
  sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
  
  # combined distances with sample data
  colnames(sd) = c("Var1", "Type1")
  wu.sd = left_join(wu.m, sd, by = "Var1")
  
  colnames(sd) = c("Var2", "Type2")
  wu.sd = left_join(wu.sd, sd, by = "Var2")
  
  #remove this line to plot all comparisons. 
  wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
  
  # plot
  ggplot(wu.sd, aes(x = Type2, y = value)) +
    theme_bw() +
    geom_point() +
    geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
    scale_color_identity() +
    facet_wrap(~ Type1, scales = "free_x") +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
    ggtitle(paste0("Distance Metric = ", m))
  
}

a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta.pdf")

Saving 7 x 7 in image

b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta.pdf")

Saving 7 x 7 in image

c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta.pdf")

Saving 7 x 7 in image

library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_allgroup_beta.pdf", width = 7, height = 10)

Taxon abundance bar plot

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",

        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",

        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

phylumGlommed <- tax_glom(ps_rare, "Phylum")

#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")

Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion

phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Taxon_barplot_t1.pdf")

Saving 7.29 x 4.51 in image

#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")

Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion

phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t2.pdf")

Saving 7.29 x 4.51 in image

#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")

Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion

phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t3.pdf")

Saving 7.29 x 4.51 in image